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John  Strain 
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University  of  California 


Abstract  The  research  supported  by  this  grant  aimed  to  solve  moving  interface 
problems  with  two  flexible  computational  modules:  semi-Lagrangian  contouring  (SLC) 
methods  which  evolved  an  interface  with  a  given  velocity  field,  and  problem-dependent 
evaluation  of  the  velocity  of  a  given  interface.  These  modules,  with  associated  engines 
for  fast  computational  geometry,  constituted  a  highly  effective  toolbox  for  evolving 
an  implicitly  represented  interface  through  arbitrary  topology  under  a  given  velocity 
functional.  The  first  module  was  extended  with  linear ly-implicit  SLC  methods  which 
evolved  stiff  interfaces  under  geometric  velocities  such  as  mean  curvature  and  surface 
diffusion,  and  fourth-order  SLC  methods  which  delivered  greatly  improved  accuracy. 
The  second  module  was  extended  with  fast  solvers  for  elliptic  systems  which  evalu¬ 
ated  interface  velocities  determined  by  quasistationary  physical  interactions  off  the 
interface,  and  enabled  the  solution  of  key  models  such  as  Ostwald  ripening,  elastic 
membranes  in  Stokes  flow,  and  crystal  growth. 


Research  results 


A  moving  interface  is  a  collection  F(f)  of  nonintersecting  oriented  closed  curves  or 
surfaces  with  an  outward  unit  normal  A,  curvature  C ,  and  normal  velocity  V,  speci¬ 
fied  as  a  functional  of  T(£).  Examples  include  passive  transport  where  V  =  F(x,t)  is 


given;  geometric  motion  V  =  ( p(N)—q(N)C)N ;  Ostwald  ripening  V  =  N  —  AC, 
where  A u  =  0  off  T(f),  u  —  C  on  T(i),  and  A  is  the  Dirichlet-Neumann  operator;  and 

N ,  where  ut  =  A u  off  T(i)  and  u  =  f(N,  C,  V )  on 
ace  problems  as  evolution  of  an  implicit  represen¬ 


ts 


dN 


models  for  crystal  growth  V  = 

T (t).  We  reformulate  moving  inter: 
tation  function  ip  whose  zero  set  is  T(t),  plus  some  densities  defined  on  T(t).  Then  at 
a  new  time  level  t  +  k,  the  new  implicit  representation  is  given  by  ip(x,t  +  h)  =  ip(x,  t) 
where  x  is  the  foot  of  the  characteristic  ending  at  x. 

Semi-Lagrangian  contouring  (SLC)  is  reviewed  in  §1.1,  applied  to  Stokes  flow  in 
§1.2,  and  to  3D  viscoelastic  flow  in  §1.3.  Two  results  in  fast  computational  geometry 
are  reported  in  §1.4,  and  employed  in  §1.5  for  fourth-order  SLC.  Linearly-implicit 
SLC  for  stiff  problems  is  presented  in  §1.6. 


1 


1.1.  Semi-Lagrangian  contouring  Semi-Lagrangian  contouring  (SLC)  delivers  a 
coherent  and  widely  usable  black-box  module,  which  moves  an  interface  through  ar¬ 
bitrary  topological  and  geometric  events,  with  minimal  information  about  the  specific 
physical  problem  and  the  interface  velocity  V.  It  incorporates  implicit  interface  evo¬ 
lution,  semi-Lagrangian  advection,  fast  geometric  algorithms  and  fast  elliptic  solvers. 

Each  point  of  a  moving  interface  T  travels  along  a  characteristic  curve  of  the 
interface  velocity  V,  carrying  a  zero  value  of  the  signed  distance 

d{x)  =  ±  min  \\x  —  7||. 

7er 

We  extend  the  velocity  V  smoothly  to  a  global  field  W,  move  all  values  of  d  along 
approximate  characteristics  of  W,  and  extract  the  resulting  interface  by  contouring 
semi-Lagrangian  approximate  advection  formulas  such  as  the  first-order  Courant- 
Isaacson-Rees  (CIR)  scheme 

ip(x)  =  d(x,  t  +  k)  ~  d{x  +  kV (x,  £),£).  (1) 

The  analogous  second-order  scheme  couples  the  first-order  predictor  (1)  with  a  trape¬ 
zoidal  corrector  based  on  the  averaged  velocity.  It  combines  the  unconditional  stabil¬ 
ity  of  the  first-order  scheme  with  the  dramatically  reduced  dissipation  of  the  trape¬ 
zoidal  rule.  Exact  distance  finding  in  a  dynamic  quadtree  data  structure  or  an  adapted 
Voronoi  diagram  eliminates  the  usual  semi-Lagrangian  interpolation  error.  The  in¬ 
terface  velocity  is  extended  by  a  numerical  Whitney  extension,  which  satisfies  a  max¬ 
imum  principle  and  removes  the  discontinuities  of  the  usual  nearest-point  extension. 
After  evolving  the  implicit  representation  ip,  the  new  interface  is  extracted  by  sub¬ 
grid  refinement  of  a  tree  mesh:  given  a  signed  distance  function  (p(x,tn),  we  build 
a  tree  at  time  tn+ 1  =  tn  +  k  by  recursive  evaluation  of  g(x)  =  <p(s,tn)  at  projected 
points  s  =  x  +  kW(x,tn).  SLC  methods  converge  to  correct  viscosity  solutions  for 
difficult  moving  interface  problems  involving  merging,  faceting,  transport,  nonlocality 
and  anisotropic  curvature-dependent  geometry  [12]. 

1.2.  Stokes  flow  with  elastic  interfaces  Slow  viscous  flows  which  satisfy  the 
incompressible  Stokes  equations 

—  vAu  +  Vp  —  F,  V  •  u  =  0 

commonly  occur  in  biological  moving  interface  problems,  for  example  blood  flow,  cell 
movement,  and  atherosclerosis.  The  force  F  in  these  equations  is  often  modeled  by  a 
measure  F  =  /hr  when  the  interface  T  has  a  complex  internal  elastic  structure.  We 
have  developed  a  fast  solver  for  Stokes  flow  with  elastic  interfaces,  by  combining  semi- 
Lagrangian  interface  evolution  with  a  fast  new  Ewald  summation  scheme  [3].  The 
semi-Lagrangian  transport  of  interface  densities  allows  a  remarkably  straightforward 
computation  of  stretching  energy,  while  the  new  fast  summation  technique  unifies 
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several  well-known  local  correction  techniques  for  singular  integral  operators.  Our 
technique  remains  highly  effective  with  discontinuous  data,  where  standard  Stokes 
solvers  encounter  substantial  difficulty.  Figure  1  exhibits  circularizing  and  oscillating 
interfaces  under  Stokes  flow,  accurate  to  one  part  per  thousand. 


Figure  1:  Circularizing  and  oscillating  interfaces  in  Stokes  flow. 


1.3.  3D  viscoelastic  flow  with  structured  interfaces  The  semi-Lagrangian 
method  has  been  implemented  in  3D  and  coupled  with  a  viscoelastic  fluid  simulator 
to  produce  complex  and  realistic  fluid  simulations  [1,  2].  It  is  publicly  available  as 
part  of  the  open-source  Berkeley  Fluid  Animation  and  Simulation  Toolkit  (BFAST) 
on  SourceForge.  The  semi-Lagrangian  approach  is  capable  of  tracking  local  features 
such  as  the  surface  colors  shown  in  Figure  2  or  the  elastic  energy  density  of  a  tensile 
membrane  [3]. 

1.4.  Computational  geometry  Operations  on  implicit  representations  are  widely 
useful  in  computational  science.  We  are  extending  our  efficient  2D  piecewise-linear 
modules  for  implicit  geometric  operations  such  as  contouring,  distancing  and  exten¬ 
sion,  to  3D  higher-order  piecewise-polynomial  interfaces. 
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Figure  2:  Paint  balls  and  merging  sprays  in  3D  viscoelastic  flow. 


Accurate  contouring  The  general  contouring  problem,  finding  a  smooth  geo¬ 
metrically  constrained  approximate  zero  set  of  a  function  which  can  be  evaluated  at 
arbitrary  points,  occurs  frequently  in  computational  science  and  requires  a  robust 
package.  An  ideal  contouring  package  would  accept  function  values  (and  derivatives 
if  available)  at  arbitrary  points,  and  produce  a  piecewise-smooth  approximation  to 
the  zero  set,  with  corners  where  necessary.  Geometric  constraints,  such  as  bounds 
on  curvature  away  from  corners,  are  vital  in  applications  like  computer-controlled 
machining.  They  pose  a  major  complication  for  existing  public-domain  contouring 
software.  We  are  finalizing  new  C/C++  packages  for  constrained  piecewise-smooth 
contouring  in  two  and  three  dimensions,  based  on  new  local  contouring  schemes  for 
Bezier  patches  and  new  methods  of  scattered  data  interpolation  [6,  7]. 

Fast  distance  to  Bezier  curves  Computation  of  distances  to  an  explicit  surface 
is  a  key  step  in  implicitization,  and  extremely  expensive  if  implemented  directly.  Fast 
distancing  algorithms  are  available  for  standard  piecewise- linear  polygonal  surfaces. 
Smoother  piecewise-polynomial  Bezier  surfaces  would  provide  higher  accuracy,  but 
computing  distances  to  these  surfaces  has  been  intractable.  The  classical  Voronoi 
diagram  provides  efficient  algorithms  for  distance  finding  to  “sites”  which  are  points  in 
space,  but  fails  when  the  sites  are  Bezier  segments  of  a  piecewise-smooth  surface.  We 
have  developed  a  fast  Voronoi-based  distance  finding  algorithm  for  Bezier  segments, 
which  dramatically  speeds  up  the  computation  of  distance  to  piecewise-polynomial 
curves  and  surfaces  [5].  Figure  3  demonstrates  the  speed  and  robustness  of  our 
algorithm,  which  requires  no  more  than  five  elementary  distance  finding  steps  per 
evaluation. 


Figure  3:  Fast  distancing  to  piecewise-polynomial  curves:  the  background  shading 
from  white  to  black  indicates  regions  where  the  algorithm  requires  1  to  5  elementary 
distancing  operations  respectively. 


1.5.  Fourth-order  accurate  SLC  methods  New  SLC  methods  based  on  piecewise- 
cubic  contouring,  fast  distancing  and  fourth-order  time  discretization  dramatically  re¬ 
duce  the  number  of  degrees  of  freedom  required  to  represent  complex  nearly-singular 
interfaces  (see  Fig.  4) [8,  10]. 


Figure  4:  Fourth-order  cubic  approximation  of  a  nearly  singular  contour  (right)  is 
far  more  efficient  than  second-order  linear  approximation  (left). 

Fourth-order  time  discretization  integrates  the  characteristic  terminal  value  prob¬ 
lem 

x'(s)  =  W(x(s),  s)  x(t  +  k)  —  x 

backward  in  time  till  s  =  t.  Input  data  consists  of  x(t  +  k)  —  x  at  the  terminal  time 
t+k,  and  W (x,  t )  at  the  initial  time  t.  Subsequent  values  of  W (x,  s )  become  available 
only  after  each  new  approximation  of  x(s)  has  been  computed.  The  mismatch  in  time 
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between  x  and  W  requires  unconventional  time  discretization  schemes  such  as  the 
first-order  CIR  predictor 


xq  =  x,  Ww  —  W(x,t),  x\  =  x  —  kW\Q 
and  the  second-order  trapezoidal  corrector 

Wio  —  W(xi,t),  Wqi  =  W(xo,t  +  k),  X2  =  x  — —W$i  — iq. 

Thus  we  derive  higher-order  correctors 

W20  =  W (x2,  t ),  x3  =  x  -  akW20  ~  (3kW10  -  ^kWm 

by  the  method  of  undetermined  coefficients  [10].  Numerical  computations  of  passive 
rotation  and  geometric  evolution  demonstrate  the  efficiency  of  fourth-order  SLC  for 
moving  complicated  smooth  interfaces  under  nonstiff  velocity  fields  (Fig.  5). 


(a)  (b)  (c) 

Figure  5:  High-order  semi-Lagrangian  contouring  evolves  smooth  and  nonsmooth 
shapes  accurately  under  smooth  geometric  velocity  fields. 


1.6.  Stiff  problems  The  interface  velocity  V  often  depends  011  the  moving  interface 
r  (t)  in  a  stiffly-stable  way.  Local  geometric  problems  such  as  mean  curvature  flow  V  = 
—CN  and  surface  diffusion  V  =  (A C)N  are  equivalent  to  stiff  nonlinear  parabolic 
PDEs  for  the  implicit  interface  representation  (p.  Nonlocal  problems  of  materials 
science  (crystallization,  Ostwald  ripening,  Hole  Shaw  flow)  determine  the  interface 
velocity  from  a  global  system  of  PDEs  with  stiff  geometric  boundary  conditions.  Our 
new  stiffly-stable  implicit  SLC  methods  and  damped  Newton-Krylov  solvers  provide 
efficient  solution  methods  for  stiff  moving  interface  problems. 

For  stiff  moving  interface  problems,  standard  explicit  formulas  such  as  first-order 
CIR  scheme 

g(x)  :=  <p(x,  t  +  k)  =  cp(x  —  kW (x,  t),t) 
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are  stable  in  the  maximum  norm  but  not  in  norms  involving  higher  derivatives.  The 
resulting  interface  oscillations  can  be  stabilized  by  taking  small  time  steps  or  frequent 
smoothing. 

Efficient  solution  of  stiff  problems  requires  implicit  time  discretizations  such  as 
the  first-order  implicit  formula 

g(x)  =  ip(x  -  EgVg(x,t  +  k),t)  =:  F[g](x) 

in  which  the  known  velocity  W (x,  t)  is  replaced  by  the  unknown  velocity  W (x,  t+k )  = 
EgVg(x,t  +  k).  The  velocity  W(x,t  +  k)  depends  on  the  unknown  implicit  interface 
representation  g(x)  ps  (p(x,t  +  k)  twice.  The  operator  Eg  extends  functions  from  the 
zero  set  T  of  g  to  the  whole  ambient  space,  while  the  problem-dependent  interfacial 
velocity  Vg  depends  stiffly  on  derivatives  of  g.  We  solve  the  implicit  interface  motion 
by  damped  Newton  iteration  and  Krylov  space  methods,  and  obtain  excellent  results 
for  mean  curvature  flow.  Nonlocal  problems  are  treated  by  the  chain  rule  and  a 
library  of  basic  Frechet  derivatives. 


Damped  Newton  iteration 

g  such  as 

N 


Since  Vg  usually  depends  on  geometric  derivatives  of 


Vg 

livin’ 


C  =  V  ■  N, 


the  simple  fixed  point  iteration  g  •<—  g  —  F\g\  will  not  converge.  Thus  we  employ 
damped  Newton  iteration 


DF[g]h  = -F[g],  g  «-  g  +  Xh. 

Here  DF[g]  is  the  Frechet  derivative  of  F  at  g.  The  damping  parameter  A  is  adaptively 
chosen  to  guarantee  steady  decrease  of  the  residual:  ||F[g]||  4—  ||F[g]||/2.  Damped 
Newton  iteration  preconditions  the  loss  of  derivatives  in  each  F  evaluation  by  the 
inverse  Jacobian  DF\g]~1,  yielding  a  modified  fixed  point  iteration  with  improved 
convergence. 

At  each  damped  Newton  step,  the  variational  equation  DF[g\h  =  —F[g\  is  approx¬ 
imated  by  a  large  sparse  system  of  linear  equations  and  solved  by  linear  or  nonlinear 
Krylov  space  methods.  Linear  methods  such  as  GMRES  converge  faster,  but  each 
step  requires  an  expensive  evaluation  of  the  Jacobian  DF .  Nonlinear  methods  require 
only  F  evaluations. 

Implicit  SLC  with  damped  Newton  iteration,  solved  by  Krylov  space  methods, 
deliver  efficient  algorithms  for  stiff  geometric  moving  interface  problems  [8,  11].  Nu¬ 
merical  results  for  mean  curvature  flows  of  smooth  and  nonsmooth  objects  are  shown 
in  Fig.  6. 


Chain  rule  The  Frechet  derivative  D^V  involved  in  second-order  elliptic  moving 
interface  problems  (such  as  Ostwald  ripening  V  =  A ^C)  is  a  problem-dependent 
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(a)  (b)  (c) 

Figure  6:  Implicit  semi-Lagrangian  contouring  evolves  both  circles  (superimposed  on 
unstable  explicit  results)  and  nonsmooth  shapes  stably  under  mean  curvature  flow. 


third-order  linear  pseudodifferential  operator.  It  can  be  constructed  from  problem- 
independent  components  by  the  chain  rule: 


DpV  —  Dlp(AlpC)  —  D(pA(p[C)  +  Atp(DipC) . 


Hadamard’s  variational  formula  implies  an  augmented  elliptic  system  for  the  first 
component  DpAp,  derived  by  implicit  differentiation  of  the  PDEs  and  boundary 
conditions  which  define  Ap.  The  second  component  DpC  is  a  problem-independent 
local  geometric  sensitivity.  Efficient  problem-independent  modules  can  compute  and 
assemble  these  components  of  the  Frechet  derivative  for  a  wide  variety  of  stiff  moving 
interface  problems. 

2.  Fast  elliptic  solvers  Elliptic  partial  differential  equations  (PDEs)  for  the  inter¬ 
face  velocity  are  converted  to  first-order  systems  in  §2.1  and  solved  by  finite  differences 
in  §2.2.  New  Alternating  Direction  Implicit  (ADI)  iterations  for  elliptic  systems  are 
presented  in  §2.3,  and  generalized  Ewald  summation  for  elliptic  systems  in  §2.4.  Two 
nonuniform  fast  Fourier  Transforms  for  piecewise-polynomial  distributions  are  de¬ 
scribed  in  §2.5  and  §2.6,  and  employed  in  §2.7  and  §2.8,  to  speed  up  locally-corrected 
boundary  integral  solution  of  elliptic  systems  in  complex  domains. 

2.1.  First-order  systems  Industrial  moving  interface  problems  (such  as  solidifica¬ 
tion,  fluid-structure  interaction  and  crystal  growth)  determine  the  interface  velocity 
V  by  solving  systems  of  PDEs  on  the  moving  phase  domains,  with  boundary  con¬ 
ditions  on  fixed  and  moving  boundaries.  Moving  interfaces  with  our  SLC  methods 
relies  on  a  problem-dependent  module  which  evaluates  the  interface  velocity  V  of 
the  interface  Y  =  {ip  =  0}.  Implicit  time  discretization  of  stiff  problems  requires 
the  Frechet  derivative  DpV  with  respect  to  the  interface  representation  ip ,  which 
satisfies  additional  PDEs.  New  locally-corrected  spectral  boundary  integral  methods 
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[13,  17,  16]  solve  general  elliptic  systems  of  PDEs  with  complex  interfaces:  Linear 
constant-coefficient  elliptic  systems  (such  as  Maxwell,  Stokes,  or  linear  elasticity)  are 
converted  to  overdetermined  first-order  systems 

d 

Au(x)  =  ^2  Ajuj(x)  +  A0u(x)  =  f(x)  in  hi,  Bu  =  g  on  T  =  dQ  (2) 
o= i 

where  each  Aj  is  a  p  x  q  matrix,  B  =  B{ 7)  is  an  r  x  q  matrix  for  each  7  6  T,  and  u 
is  a  q- vector.  Such  systems  are  amenable  to  new  iterations  and  simple  new  boundary 
integral  solvers. 

2.2.  Finite  difference  methods  A  new  piecewise-polynomial  interface  method 
for  discretizing  elliptic  systems  with  complex  interfaces  between  high-contrast  mate¬ 
rials  is  derived  and  analyzed  in  [4] .  A  Krylov-accelerated  interface  multigrid  approach 
solves  the  new  discretization  efficiently.  Stability  and  convergence  are  proved  in  one 
dimension.  Numerical  experiments  with  complex  two-dimensional  interfaces  (see  Fig. 
7)  and  coefficients  varying  over  eight  orders  of  magnitude  demonstrate  the  accuracy, 
efficiency  and  robustness  of  the  method. 


y  -1  -1  x  y  -1  -1  x 

(c)  (d) 

Figure  7:  Solution  of  elliptic  problems  with  complex  interfaces. 
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2.3.  ADI  methods  for  elliptic  systems  Classical  ADI  methods  provide  the 
oldest  essentially  optimal  iteration,  for  separable  second-order  positive  definite  elliptic 
problems 

Au  +  Bu  =  —dx  ( adxu )  —  dy  ((3dyu)  =  /, 

in  rectangular  two-dimensional  domains  C  R2.  With  scale  1/s,  these  iterations 
take  the  form 

(s2  +  A)  (s2  +  B )  um+1  =  (s2  -  A)  (s2  -  B)  um  +  2/ 

and  modify  each  Fourier  mode  elkTx  with  wavenumber  k  =  (ki,k2)  by  a  symbol  a 
which  damps  error  rapidly  over  a  geometric  range  of  wavenumbers.  For  example,  in 
the  Poisson  equation  |er(/c)|  <  1  whenever  -j=  <  ^  <  y/2.  Thus  varying  the  scale  s 
geometrically  gives  0(e)  damping  in  0(log  Aloge)  sweeps.  Each  sweep  approximates 
(A  +  B)^1  by  a  product  of  one- dimensional  operators  (s2  ±  A)±1  and  requires  time 
proportional  to  the  problem  size,  making  ADI  an  essentially  optimal  iteration. 

We  derive  ADI  methods  for  general  elliptic  systems  (2).  Ellipticity  implies  that 
in  any  bounded  domain  with  boundary  Y  =  dfl  and  inward  normal  vector  field  n, 
the  inward  normal  derivative  of  u  is  determined  by  values  and  tangential  derivatives 
of  u  on  T: 

dnu  =  Y  ni,9iu  =  Al(f  -  ATdTu  -  A0u) .  (3) 

i 

Here  the  left  inverse  A\  of  An  =  Xq=1  rijAj  exists  by  ellipticity,  and  At&t  is  the  tan¬ 
gential  part  of  the  operator  A  =  Andn  +  AtBt ■  Given  the  inward  normal  derivative, 
the  solution  can  be  approximated  by  Taylor  expansion  on  a  strip  near  the  boundary. 
Marching  inward  with  repeated  Taylor  expansion  solves  the  boundary  value  problem, 
and  explains  heuristically  why  boundary  value  problems  are  well-posed  for  elliptic 
systems.  Thus  our  ADI  methods  sweep  inward  from  the  boundary  to  provide  conver¬ 
gent  three-step  iterations  for  elliptic  systems:  1.  Choose  sweep  direction  n  pointing 
in  from  the  boundary  and  rewrite  the  system  in  normal/tangential  coordinates  (3). 

2.  Left-invert  An  and  shift  by  the  inverse  length  scale  s  to  get  the  equivalent  system 

su  +  dnu  +  Bq  u  =  su  +  f  —  Bt8tu  +  BqU.  (4) 

3.  Solve  Eq.  (4)  inward  from  the  boundary  across  the  domain,  with  step  size  pro¬ 
portional  to  1/s.  Repeating  these  three  steps  with  a  sequence  of  directions  n  and 
geometrically  varying  scales  s  gives  a  full  ADI  iteration. 

Convergence  In  rectangular  domains  with  periodic  boundary  conditions,  a  col¬ 
lection  of  sweeps  with  different  directions  n  reduces  each  error  mode  QlkTxJ  by  a 
matrix-valued  symbol 

cr(k)  =  IJ  (s  +  lkn  +  BoY  (s  ~  iBTkT) . 
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Employing  a  larger  set  of  directions  improves  the  convergence  speed  exponentially: 
one  pass  of  ADI  iteration  reduces  the  error  by  0.5  with  4  directions,  0.1  with  8 
directions,  and  0.005  with  16  directions.  Since  this  error  reduction  is  independent  of 
mesh  size  (see  Figure  8),  the  new  ADI  iteration  is  an  essentially  optimal  solver  for 
elliptic  systems  in  simple  domains  [14]. 


4  directions 


8  directions 


16  directions 


Figure  8:  ADI  error  reduction  symbol  a  versus  wavenumber  k  for  the  Poisson  equa¬ 
tion  on  a  IV  =  2562  mesh,  solved  with  4,  8  or  16  directions  in  0(N  log N)  work. 


2.4.  Generalized  Ewald  summation  For  elliptic  systems  with  complex  inter¬ 
faces,  locally-corrected  spectral  methods  [13,  17,  16]  employ  the  periodic  fundamental 
solution 

S{x-y)  =  J2  (e  iM.  +  Aj)  eikT(x~y)  =  J2  a{k)Wr(x-v\  af  -  (a*a)"V 

k£Zd  \j= 1  /  kezd 

of  the  first-order  system  (2).  Generalized  Ewald  summation  [13]  extracts  a  global 
rapidly-converging  Fourier  series 

ST(x-y)=  J2  e_Ta*(fc)a(fc)a(/c)teifcT^_y\ 
kezd 

mollified  by  a  matrix  exponential  e~Ta*  (fcMfc) .  The  error  term  has  a  formal  asymptotic 
expansion  as  r  — >  0 

r2 

Er  =  S  -  ST  =  (I  -  e~TA*A)(A*A)~1A *  =  tA*  -  —(A*A)A*  +  •  •  •  (5) 

which  corrects  via  an  asymptotic  series  of  local  differential  operators.  These  new 
mollification  and  local  correction  techniques  combine  with  the  fast  Fourier  transform, 
Pade  codes  for  small  dense  matrix  exponentials,  and  high-order  uncentered  differ¬ 
encing  to  solve  first-order  elliptic  systems  in  periodic  geometry.  A  simple  algebraic 
algorithm  automatically  computes  local  correction  coefficients  which  achieve  high- 
order  accuracy  at  minimal  cost. 
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Locally-corrected  spectral  volume  potentials  Fast  solvers  for  elliptic  systems 
(2)  employ  the  volume  potential 

Vf(x)  =  [  S(x  —  y)f(y)  dy 
Jci 

to  eliminate  the  right-hand  side  /.  Such  potentials  solve  elliptic  systems  with  piecewise- 
smooth  right-hand  sides  which  are  discontinuous  across  the  interface  T  =  dfl.  Locally- 
corrected  spectral  methods  treat  the  interfacial  discontinuities  as  distributions,  con¬ 
verting  the  problem  to  an  inhomogeneous  system  Au  =  F  with  a  singular  right-hand 
side  F  =  fxn  on  a  fixed  enclosing  domain  B  D  Q.  Here  /  is  smooth,  while  the 
characteristic  function  °f  the  domain  jumps  from  1  to  0  across  the  interface 
T  and  includes  all  the  singularities  of  the  data.  Generalized  Ewald  summation  lo¬ 
calizes  and  separates  the  singular  interface  distributions  to  solve  the  inhomogeneous 
system  accurately  and  efficiently.  The  Fourier  series  VT  =  ST  *  F  is  computed  by  a 
geometric  nonuniform  fast  Fourier  transform  (GNUFFT)  [9,  15].  The  local  correction 
(5)  separates  smooth  from  nonsmooth  parts  of  the  solution  by  the  product  rule  for 
derivatives.  For  example,  the  first-order  correction  Er  =  tA*  +  0(r2)  yields 

rA*(fxn)  =  t  ^J2A*(fjXn  +  frijSr)  +  A *0fxnj 

where  n  is  the  outward  unit  normal  to  fl  and  is  a  delta- function  on  F.  Since  the 
enclosing  domain  is  fixed  and  simple,  the  resulting  locally-corrected  volume  potential 
capitalizes  on  fast  solvers  for  simple  problems  to  achieve  optimal  efficiency  [17]. 


2.5.  GNUFFT  by  B-spline  smoothing  Fast  solvers  based  on  Ewald  summa¬ 
tion  require  the  computation  of  accurate  low-frequency  Fourier  coefficients  of  sin¬ 
gular  distributions  such  as  5-functions  spread  over  curves,  surfaces,  and  other  lower¬ 
dimensional  geometric  objects.  Standard  nonuniform  FFT  (NUFFT)  algorithms  work 
only  for  point  distributions.  New  geometric  NUFFTs  or  GNUFFTs  compute  accurate 
Fourier  coefficients  of  a  piecewise-polynomial  distribution  supported  on  arbitrary- 
dimension  simplices  scattered  in  D- dimensional  Euclidean  space.  Each  singular  term 
is  mollified  with  a  B-spline  smoothing  kernel,  evaluated  on  a  nearby  uniform  mesh, 
transformed  with  a  standard  FFT,  and  deconvolved  in  real  space.  The  resulting 
algorithm  displays  dramatic  speedups  over  direct  evaluation,  reduces  to  the  stan¬ 
dard  NUFFT  in  simplex  dimension  zero,  and  provides  a  highly  effective  tool  for  the 
locally-corrected  spectral  solution  of  elliptic  systems  [9]. 


2.6.  GNUFFT  via  the  butterfly  algorithm  A  new  GNUFFT  based  on  the 
butterfly  algorithm  generalizes  pointwise  NUFFTs  based  on  low-rank  approximations 
such  as  the  Taylor  expansion 


etiSj  =  etia  V  ^ — ^e~TaeTSi  +  0(pm/m\), 
~‘n  al 

a= U 
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valid  in  any  Heisenberg  rectangle  1 1  —  r||s  —  a\  <  p.  The  butterfly  GNUFFT  builds 
moments  due  to  sources  Sj  in  Heisenberg  rectangles  with  minimal  span  in  source 
space,  recursively  merges  and  shifts  them  to  maximize  source  span  and  minimize 
target  span,  and  evaluates  the  final  expansions  at  targets  in  optimal  0(N  log  N)  time 
[15].  The  recursion  is  illustrated  in  Fig.  9. 


Figure  9:  Recursive  subdivision  of  a  spiral  of  simplices  for  the  butterfly  GNUFFT. 


2.7.  Locally-corrected  spectral  boundary  integral  solvers  A  boundary  in¬ 
tegral  equation  is  derived  in  [13]  with  the  periodic  fundamental  solution  S(x  —  y)  of 
the  elliptic  system  (2):  multiplying  Eq.  (2)  by  S,  integrating  over  H  and  applying 
Gauss’  theorem  yields 

[  Pil) S(l-cr) An{a)n(a)  da  =  p(~/),  An(a)  =  Yjnj(a)Aj.  (6) 

z  Jr  j= i 

The  new  unknown  =  P(p/)u(^)  =  (/  —  B*(^)B{^))u{^)  is  the  projection  of  u 
orthogonal  to  the  data  g.  The  right-hand  side  is  a  combination  of  volume  and  layer 
potentials 

P{l)  =  [  ~y)f(y)dy+  I  -  a)An(a)B* g(a)  da. 

Jn  Jr 
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Eq.  (6)  is  square,  well-conditioned,  and  amenable  to  solution  by  a  variety  of  fast  stable 
high-order  accurate  methods.  The  solution  u  of  the  elliptic  system  (2)  is  recovered 
locally  on  the  boundary  by  u( 7)  =  /i( 7)  +  B*(nf)g('y),  and  globally  by  integration. 

The  boundary  integral  equation  (6)  is  solved  with  generalized  Ewald  summation 
for  S  (§2.4),  which  approximates  the  kernel  if  (7,  cr)  =  P(7)S'(7  —  a)An(a)  by  a 
global  rapidly-converging  Fourier  series,  and  converts  the  integral  equation  into  semi- 
separated  form 

(/  -  MRT)p  =  p.  (7) 

Here  Tp  evaluates  Fourier  coefficients  of  (Hnp)hr,  R  is  the  hltered  inverse  of  the 
elliptic  operator  in  Fourier  space,  and  M  evaluates  and  projects  Fourier  series  on  T. 
Since  (/  —  MRT)~l  =  I  +  MR(I  —  T M R)~1T ,  the  semi-separated  form  (7)  can  be 
solved  in  Fourier  space,  where  the  matrix  of  TM R,  is 

f  (k,q)  =  2  J^An{a)e~'l{k-q)Ta  dae~Ta^q)a^ a* (q) . 

This  matrix  contains  the  elliptic  system,  interface  and  boundary  conditions.  Its 
regular  data  structure  is  amenable  to  fast  randomized  low-rank  approximation  and 
efficient  solution.  Local  corrections  are  higher-order  in  the  small  parameter  r,  and  can 
be  applied  locally  after  each  Fourier  space  iteration.  The  resulting  boundary  integral 
solver  is  efficient,  stable,  and  as  accurate  as  the  underlying  interface  representation. 

Frechet  derivatives  Our  boundary  integral  approach  expresses  the  interface  ve¬ 
locity  V  by  a  sequence  of  three  explicit  differentiable  steps  which  facilitate  the  ap¬ 
proximation  of  Frechet  derivatives.  First,  the  solution  u  of  the  elliptic  system  is  a 
local  function 

M(7)  =  (I  -  B*B)u( 7)  +  B*g( 7)  =  /i(y)  +  B*g( 7) 

of  /i  on  the  interface.  Its  values  off  the  interface  depend  only  on  the  known  values  on 
the  interface.  Second,  the  integral  density  p  satisfies  the  boundary  integral  equation 

(/  -  MRT)p  =  p, 

where  the  operators  M  and  T  depend  on  ip.  The  damped  pseudoinverse  R  of  the 
elliptic  operator  A  is  independent  of  the  interface.  Third,  the  right-hand  side  p  of 
the  boundary  integral  equation  depends  on  p  via  integration  over  the  interface  and 
projection  by  the  boundary  matrix  B.  In  moving  interface  problems,  g  depends  on  the 
normal  n  and  curvature  C ,  and  DpV  involves  the  corresponding  Frechet  derivatives. 

Given  this  three-step  construction  of  the  solution  u,  the  standard  formula  D(X^1)  = 
—X~1DXX~1  and  the  chain  rule  exhibit  Dpp  as  a  sum 

Dyii  =  (/  -  MRT)-1  ((D^M)RT  +  MRD^t)  (/  -  MRT)~lp  +  (/  -  MRT)~1Dlpp 

(8) 

of  products  of  complicated  operators.  Each  factor  encapsulates  the  linearized  response 
of  one  component  of  the  system  to  a  perturbation  in  the  interface.  Eq.  (8)  determines 
the  exact  Frechet  derivative  D^pV  and  generalizes  the  classical  Hadamard  formula. 
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3.  Interactions/Transitions 

The  PI  presented  results  from  this  research  in  seminars  and  colloquia  at  FAN2010: 
Fluid  dynamics,  analysis  and  numerics  (Duke  University),  the  IMA  Hot  Topics  Work¬ 
shop  on  Integral  Equation  Methods,  Fast  Algorithms  and  Applications  (Minneapo¬ 
lis),  the  International  Conference  on  Industrial  and  Applied  Mathematics  (Zurich 
and  Vancouver  ),  the  International  Conference  on  Advances  in  Scientific  Computing 
(Brown  University),  California  Institute  of  Technology,  Duke  University,  Michigan 
State  University,  the  National  University  of  Singapore,  North  Carolina  State  Uni¬ 
versity,  Purdue  University,  University  of  North  Carolina  at  Charlotte,  University 
of  California  at  Berkeley,  University  of  California  at  Irvine,  the  American  Institute 
of  Mathematics,  the  Institute  for  Mathematics  and  its  Applications,  INRIA  Paris- 
Rocquencourt,  and  the  Statistical  and  Applied  Mathematical  Sciences  Institute. 

Efficient  and  accurate  2D  contouring  and  distancing  codes  are  nearing  comple¬ 
tion,  and  will  soon  be  open-sourced  from  the  Pi’s  website  [5,  7].  High-order  implicit 
semi-Lagrangian  contouring  codes  for  moving  interfaces  will  follow  shortly  thereafter 
[10,  11].  Geometric  nonuniform  fast  Fourier  transform  codes  [9,  15]  and  fast  boundary 
integral  solvers  for  elliptic  systems  [17,  16,  14]  are  in  progress.  The  3D  surface  track¬ 
ing  method  described  in  [1,  2]  has  been  implemented  as  part  of  the  Berkeley  Fluid 
Animation  and  Simulation  Toolkit  (BFAST),  which  has  been  open  source  released 
and  is  available  on  SourceForge. 
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